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Abstract 

Quantum  mechanical  calculations  based  on  density 
functional  theory  (DFT)  are  used  to  study  dynamic 
behavior  of  shocked  polymeric  nitrogen,  a novel  energetic 
material.  We  report  results  on  system  sizes  in  excess  of 

3.000  atoms.  Such  calculations  on  system  sizes  within  the 

1.000  atom  range  remain  problematic  using  standard 
implementations  of  DFT.  We  evaluate  the  feasibility  of 
using  several  available  DFT  codes  for  this  work  through 
comparison  of  scalability  and  resource  requirements.  In 
this  study,  we  utilize  a recently  developed  highly-scalable 
localized  orbital  DFT  code,  CP2K,  designed  to  treat  large 
systems.  Scaling  and  performance  benchmarks  of  the 
CP2K  on  several  Department  of  Defense  (DoD)  high 
performance  computing  (HPC)  computers  are  presented 
for  a variety  of  system  sizes  and  shapes.  Additionally,  we 
report  preliminary  calculations  on  the  conventional 
explosive  nitromethane.  In  those  calculations  in  excess  of 
3,500,  atoms  are  treated. 

1.  Introduction 

The  changing  nature  of  international  conflicts 
requires  the  transformation  of  war-fighting  capabilities 
within  the  Department  of  Defense.  The  attainment  of  the 
transformation  of  the  military  force  structure  depends  on 
the  development  of  lighter-weight,  more  robust  and 
integrated  armor  materials  and  for  systems  to  meet  the 
requirements  for  speed,  mobility,  and  rapid  deployment  of 
weapon  platforms  and  to  enhance  combat  personnel 
effectiveness,  and  survivability  in  the  face  of  a widening 
range  of  more  lethal  and  countering  threats.  New 
energetic  materials,  with  substantially  enhanced 
performance  and  reduced  vulnerabilities,  must  be 
developed  for  future  weapons  systems  to  meet  these 
criteria.  Advanced  Energetic  Materials  (AEM)  are 
required  to  enable  high  priority  military  missions  ranging 


from  Hard  and  Deeply  Buried  Target  Defeat,  to  Advanced 
Propulsion,  to  lightened  highly  mobile  force  evolution  and 
the  thrust  towards  miniaturized  munitions  and  systems.  It 
is  recognized  that  weapons  superiority  is  dependent  on  the 
development  of  AEM.  The  National  Advanced 
Energetics  Initiative  recognizes  that  developments  in 
computational  chemistry  and  physics-based  modeling 
using  HPC  are  among  areas  that  can  provide  the  key 
factors  that  will  enable  breakthroughs  in  the  performance 
of  energetic  materials.  Our  continued  efforts  focus  on  a 
fundamental  understanding  of  the  dynamic  behavior  of 
new  and  existing  energetic  materials  (EMs)  that  can  be 
used  to  design  new  or  modify  existing  weaponry  for 
optimal  performance.  In  particular,  we  are  interested  in 
obtaining  a first-principles  description  of  the  atomic-level 
details  of  shock  initiation  of  AEMs.  In  prior  years, 
atomic-level  simulations  of  shocked  materials  have  been 
limited  to  classical  molecular  dynamics  (MD)  treatments. 
While  such  simulations  have  been  extremely  useful  in 
providing  atomic-level  detail  into  dynamic  response  of 
shocked  materials,  the  classical  models  are  not  accurate  or 
physically  realistic  over  the  wide  range  of  thermodynamic 
conditions  present  under  shock  compression.  On  the 
other  hand,  the  quantum  mechanical  approach  of  Density 
Functional  Theory  (DFT)1'1  has  demonstrated  accuracy 
over  a wide  range  of  temperatures  and  pressures  and  does 
not  suffer  the  limitations  of  classical  models121.  For  this 
reason,  we  have  chosen  to  apply  DFT  MD  methods  to  the 
study  of  AEM  to  obtain  the  desired  first-principles 
description  of  material  response  to  shock  initiation. 

In  our  applications  of  DFT  to  AEM,  we  first  apply 
these  methods  on  a simple  system  as  a proof  of  principle. 
We  have  chosen  for  our  test  system  a polymorphic  phase 
of  nitrogen  called  cubic  gauche  nitrogen  (cg-N)[3l  While 
nitrogen  seems  far-removed  from  conventional  EMs 
(which  are  typically  large,  polyatomic  organic  molecular 
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crystals),  under  pressure  it  has  been  shown  to  take  on  a 
polymeric  form.  The  cg-N  form  is  a bulk  three 
dimensional  covalent  network  like  diamond,  in  which 
every  atom  has  three  neighbors  compared  to  diatomic 
molecular  forms  of  nitrogen.  Since  diatomic  nitrogen  has 
the  strongest  bond  in  nature,  the  amount  of  energy 
released  from  the  polymeric  form  as  it  transitions  to  the 
more  stable  standard  state  is  extremely  large.  The 
products  of  the  conversion  are  very  simple;  the  solid, 
which  is  composed  of  singly-bonded  nitrogen,  transitions 
to  the  diatomic  triply-bonded  gaseous  form  and  small 
chains.  The  amount  of  energy  available  for  release 
through  complete  conversion  is  several  times  that  of  RDX 
by  weight  and  volume.  Simulation  has  confirmed  that  cg- 
N is  stable  at  low  temperatures  and  pressures,  but  will 
react  very  rapidly  and  release  large  quantities  of  energy 
when  sufficiently  perturbed14'.  The  simplicity  and 
reactivity  of  this  single-element  system  make  it  an  ideal 
physically  realistic  test  case  for  simulation  of  shock 
initiation. 

In  this  'paper,  we  will  describe  our  initial  efforts  to 
study  the  response  of  shocked  cg-N  using  DFT.  This 
paper  will  briefly  describe  the  theoretical  methods,  the 
computational  details,  software  explored  to  enable 
adequate  treatment  of  the  calculations,  HPC  resources 
used,  and  results  of  the  calculations  to  date. 

2.  Level  of  Theory 

DFT  is  not  considered  the  most  accurate  quantum 
mechanics  method  available,  but  it  has  become  a state-of- 
practice  in  the  scientific  community  due  to  its 
computational  efficiency  relative  to  other,  more  accurate 
methods.  For  example,  the  CCSD(T)'5'  quantum 
chemistry  techniques  are  formally  0(N7)  where  N reflects 
the  number  of  electrons  in  the  system  in  sharp  contrast  to 
DFT  methods  which  typically  range  from  O(N)  - 0(N3). 
In  the  DFT  formalism  each  electron  is  treated 
independently  and  interacts  with  other  electrons  not 
directly  but  only  through  an  effective  potential  which  is  a 
function  of  the  total  charge  density.  This  mean-field 
approximation  drastically  reduces  the  computational  effort 
compared  to  higher  level  quantum  chemistry  methods 
while  maintaining  accuracy  sufficient  for  many 
applications.  Advances  in  approximations,  efficient 
mathematical  algorithms,  and  implementations  on 
scaleable  platforms  have  further  increased  application  of 
DFT  to  large  and  complex  systems  of  interest  to  the  DoD. 

In  this  study,  we  have  utilized  the  Perdew-Burke- 
Emzerhof'61  form  of  the  Generalized  Gradient 
Approximation  of  the  DFT  as  implemented  in  the  local 
orbital  basis  code  CP2K.  We  have  used  the  double-zeta 
valence  polarization  basis  which  corresponds  to  a total  of 
thirteen  basis  functions  per  nitrogen  atom.  Total  energies 


and  forces  were  converged  to  1.6e-4  and  2.6e-5  a.u., 
respectively  and  were  used  throughout  these  calculations. 

3.  Computational  Details 

Two  sets  of  simulations  were  performed.  The  first  is 
a DFT  MD  simulation  of  a filament  of  material  subjected 
to  flyer-plate  impact.  This  is  a non-equilibrium  molecular 
dynamics  (NEMD)  simulation  in  which  the  dynamic 
response  of  the  material  can  be  directly  observed  through 
monitoring  the  time  progression  of  atomic  positions  and 
velocities.  The  second  set  also  utilizes  DFT-MD,  except 
the  simulations  calculate  properties  of  the  material  under 
thermodynamic  equilibrium  at  various  conditions. 
Equation-of-State  (EOS)  information  is  obtained  from 
these  simulations  and  used  to  generate  the  shock 
Hugoniot,  a special  set  of  points  satisfying  the 
conservation  equations  of  mass,  momentum  and  energy 
across  the  shock  front  between  the  quiescent  crystal  and 
its  final  shocked  state.  We  will  refer  to  this  method  as  the 
Erpenbeck'7'  approach,  named  after  an  author  who  used 
classical  MD  to  generate  the  shock  Hugoniot  of  a simple 
model  material. 

A.  Direction  Shock  Simulation 

The  simulation  cell  consists  of  a filament  of  energetic 
material  composed  of  2,304  atoms  arranged  in  the 
equilibrium  cg-N  configuration  at  T=250K,  P=1  atm. 
Periodic  boundary  conditions  are  imposed  in  the  two 
directions  perpendicular  to  the  direction  of  shock 
propagation.  In  this  study,  the  shock  wave  will  initiate  at 
one  edge  of  the  filament  using  a small  segment  of  material 
in  which  each  atom  is  assigned  large  initial  velocities  in 
the  direction  of  shock  impact.  The  remaining  atoms  in  the 
filament  are  at  a thermodynamic  equilibrium  (T=250K, 
P=1  atm).  The  resulting  integration  of  the  microcanonical 
equations  of  motion  using  DFT-MD  for  this  set  of  initial 
conditions  will  simulate  shock  initiation  through  flyer- 
plate  impact.  The  position  of  the  shock  discontinuity 
through  the  filament  is  monitored,  and  when  it  approaches 
the  end  of  the  filament,  additional  material  is  added  at  the 
end  of  the  filament  ahead  of  the  progressing  shock  wave. 
The  state  of  the  material  that  is  added  is  also  in 
thermodynamic  equilibrium  and  its  atomic  arrangements 
are  consistent  with  that  of  cg-N  at  that  thermodynamic 
state. 

Since  a shock  wave  travels  faster  than  the  sound 
speed  of  the  material,  there  is  no  need  to  simulate  large 
portions  of  undisturbed  material  ahead  of  the  shock  front; 
rather,  we  have  found  in  earlier  classical  MD  studies  that 
adding  material  ahead  of  the  propagating  shock  front  is 
sufficient  to  adequately  describe  the  dynamics  of  the 
process.  This  method  has  been  used  in  many  classical 
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MD  simulations  of  shock  waves181.  Obviously,  in  such  a 
simulation  scheme,  the  simulation  size  grows  as  the  shock 
wave  progresses  through  the  material,  significantly 
increasing  the  computational  demands  during  the 
simulation. 

B.  Erpenbeck  Method'71 

The  thermodynamic  quantities  of  a material  in  the 
quiescent  and  final  shocked  states  are  related  by  the 
conservation  equations  of  mass,  momentum,  and  energy 
across  the  shock  front  and  are  represented  by  the  shock 
Hugoniot,  a set  of  EOS  points  that  satisfy  the  Hugoniot 
function: 

Ug(T,V)  = 0 = E-Eo-'A(P  + Po)(Vo-V)  (1) 

In  this  equation,  E,  P,  and  V are  the  internal 
energy/unit  mass,  pressure,  and  volume/unit  mass, 
respectively.  E0,  Po,  and  V0  are  the  internal  energy/unit 
mass,  pressure  and  volume/unit  mass  in  the  quiescent 
state.  The  procedure  for  identifying  points  on  the  shock 
Hugoniot  is  straightforward.  For  each  volume,  a series  of 
DFT-MD  simulations  of  the  Number  of  particles,  Volume, 


and  Temperature  (NVT)  ensemble  are  performed  over  a 
range  of  temperatures.  Once  convergence  of  the 
thermodynamic  properties  is  reached  for  each  simulation, 
the  Hugoniot  function  is  evaluated.  The  resulting  set  of 
calculations  will  produce  a series  of  Hg(7V,  which  is  then 
fitted  to  a polynomial  from  which  the  Hugoniot 
temperature  (THg)  can  be  identified  from  Hg(7’Hgfr  =0. 
Once  this  has  been  done,  the  dependence  of  either  the 
pressure  or  the  internal  energy  on  temperature  can  be 
determined  through  fitting  these  quantities  to  polynomials 
and  evaluating  them  at  ^Hg-  Each  of  the  calculations  used 
to  construct  the  shock  Hugoniot  was  composed  of  256 
atoms,  with  their  initial  arrangement  consistent  with  that 
of  the  cg-N  structure  at  T=298K,  P=1  atm.  Each 
trajectory  was  first  equilibrated  for  4,000  time-steps  and 
thermodynamic  values  were  averaged  over  the  subsequent 
4,000  time-steps  (1  time-step=l  fs). 

4.  Suitability  of  Available  Codes 

For  this  Challenge  Project,  we  evaluated  six  packages 
that  implement  DFT  for  solid  state  applications  (Table  1). 


Table  1.  Evaluated  DFT  Packages'9' 
Maximum 


Package 

ABINIT 

Representation 

PW 

Size 

(Atoms) 

100 

Scalability 

(Processors) 

8 

License 

Open  Source 

CP2K 

LO 

4,000+ 

512 

Open  Source 

Espresso 

PW 

1,000 

512 

Open  Source 

PARATEC 

PW 

1,000 

512 

Limited 

SIESTA 

LO 

1 ,000+ 

128 

Commercial 

Socorro 

PW 

1,000 

512 

Open  Source 

VASP 

PW 

500 

128 

Commercial 

These  differ  mainly  in  how  they  represent  the 
electrons  and  can  be  categorized  into  two  main  families: 
Planewaves  (PW)  and  Localized  Orbitals  (LO).  The  PW 
method  is  a simple  robust  formalism  with  a long  history. 
It  is  the  “gold  standard”  in  DFT  calculations  owing  to  its 
systematically  improvable  basis  (sometimes  referred  to  as 
PW  accuracy).  Its  main  drawback  is  the  large  number  of 
PWs  required  to  adequately  describe  the  electronic  space. 
Consequently,  the  number  of  atoms  that  can  be  studied 
with  available  implementations  at  this  time  is  limited  to 
~ 1,000  due  to  the  CPU  and  memory  requirements 
associated  with  large  numbers  of  PWs.  The  LO  method 
generally  offers  excellent  accuracy  at  a fraction  of  the 
computational  cost  of  the  PW  approach,  although  it  does 


not  have  a systematically  improvable  basis.  Table  1 
provides  a description  of  the  scalability  and  applicability 
of  various  DFT  codes  on  large  systems  of  atoms  and  their 
availability.  In  our  work,  we  are  primarily  interested  in 
monitoring  the  temporal  behavior  for  systems  in  the 
1,000+  atom  range.  Therefore,  PWs  do  not  appear  to  be 
an  option  for  use  at  this  time  for  these  large  system  sizes. 
Additionally,  our  primary  calculations  are  NEMD  of 
extreme  energy-releasing  events,  and  thus  PW  accuracy  is 
not  essential  as  errors  in  the  trajectory  integration  will  be 
on  the  same  order  as  those  of  LOs.  Therefore,  the  more 
appropriate  DFT  approach  use  in  this  work  is  the  LO 
method.  The  results  shown  in  Table  1 indicate  that  for 
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efficiency,  scalability,  and  system  size,  CP2K  is  the  most 
suitable  software  package  for  our  use. 

5.  Results 

A.  DFT  MD  Simulations 

Two  types  of  MD  calculations  were  performed  to 
provide  insight  into  the  nature  of  shocked  cg-N.  The  less 
computationally  demanding  is  the  generation  of  the  shock 
Hugoniot  using  equilibrium  MD  methods.  This 
methodology  does  not  require  large  numbers  of  atoms  and 
uses  MD  to  generate  the  EOS  for  the  materials.  Its 
computational  cost  comes  from  the  number  of  EOS  points 
needed  to  map  out  the  shock  Hugoniot.  Each  point  on  the 
shock  Hugoniot  (Figure  1)  was  obtained  by  interpolating 
the  Hugoniot  function  for  four  EOS  points.  In  total  this 
has  required  in  excess  of  500,000  CPU  hours  on  the 
hpcll  (SGI  Origin  3900). 


Figure  1.  Principal  shock  Hugoniot  of  cg-N  in 
pressure-density  (left)  and  pressure-temperature 
(right)  planes 


As  evident  in  Figure  1,  no  anomalous  behavior  was 
found  in  the  principal  shock  Hugoniot.  However, 
additional  NVT  simulations  showed  that  cg-N  will 
undergo  a phase  transition  at  approximately  P = 50  GPa, 
T = 4,000K;  the  material  transitioned  from  the  polymeric 
solid  form  to  liquid  composed  of  both  atoms  and 
molecules.  The  material  in  this  state  is  not  well 
characterized.  Figure  2 shows  three  snapshots  of  the 
system  at  various  time  steps  during  the  trajectory 
integration.  The  left-most  frame  shows  the  system  in  the 
initial  cg-N  arrangement,  and  the  right-most  frame  shows 
a distribution  of  atoms  and  clusters  ranging  in  size  from 
two  to  ten  atoms.  The  structure  of  the  material  is 
consistent  with  a liquid.  The  dynamic  behavior  of  the 
material  during  this  portion  of  the  trajectory  shows  that 
the  clusters  are  not  stable  as  the  atoms  rapidly  change 
partners.  It  is  not  known  whether  any  true  cluster 
formation  is  obtained.  The  middle  frame  is  representative 
of  a transient  state  of  the  material  that  appears  to  have 
some  structure.  Further  calculations  are  needed  in  order 
to  make  a definitive  characterization  of  the  transient  and 
final  states  of  the  material. 


The  more  direct  approach  to  characterizing  a shocked 
material  is  a MD  simulation  of  a shock  wave  propagating 
through  a material.  These  calculations  require  orders  of 
magnitude  more  computational  resources  than  for  those  of 
the  Erpenbeck-type  calculations  of  the  shock  Hugoniot 
described  heretofore  due  to  requirements  of  a larger 
number  of  atoms  and  the  rapidly  changing  electronic 
environment  traveling  with  the  shock  wave.  The  upper 
pane  of  Figure  3 shows  the  left  half  initial  filament  used  in 
the  shock  simulation.  The  filament  is  composed  of  2,304 
atoms  in  the  equilibrium  cg-N  arrangement.  The  shock 
wave  is  initiated  at  the  far  left  edge  by  assigning  the  432 
left-most  atoms  in  the  filament  with  an  initial  velocity  in 
the  direction  of  the  shock  propagation  corresponding  to  1 0 
km/sec.  The  108  right-most  atoms  in  the  filament  are 
frozen  during  the  dynamics  simulation  to  preclude  surface 
effects.  As  the  shock  wave  progresses,  quiescent  material 
will  be  inserted  immediately  to  the  left  of  the  frozen 
portion  of  the  filament.  In  monitoring  the  dynamic 
behavior  of  the  shocked  material,  we  find  that  it  appears 
to  be  a stiff  material.  By  this,  we  mean  that  upon  shock 
wave  passage,  the  material  does  not  maintain  significant 
deformation  and  reverts  to  its  near  equilibrium 
configuration.  We  note  that  well  after  shock  wave 
passage,  the  atoms  at  the  left-most  edge  of  the  filament 
exhibit  large  vibrational  amplitude  motion  and  the 
material  at  this  edge  begins  to  decompose.  The  trajectory 
has  not  yet  progressed  to  the  point  that  we  observe  final 
product  formation;  however,  formation  of  the  molecular 
species  appears  to  occur,  as  shown  in  the  lower  frame  of 
Figure  3,  which  gives  a snapshot  of  the  shocked  material 
1.0  ps  after  flyer  plate  impact.  The  material  that  has 
decomposed  from  the  filament  appears  to  exist  as  either 
single  atoms  or  hot  di-  or  triatomics.  Up  to  this  point,  the 
calculations  of  the  direct  shocked  filaments  using  128 
processors  on  the  John  von  Newmann  cluster  (JVN)  have 
expended  1,000,000  CPU  hours.  Additionally  200,000 
CPU  hours  were  expended  on  hpcll  with  up  to  256 
processors. 

B.  Scaling  and  Performance  of  CP2K 

In  this  portion  of  the  paper,  we  will  describe  the 
scaling  and  performance  of  the  DFT  code  CP2K.  “Fixed” 
denotes  timings  for  a same-size  problem  across  N 
processors.  “Scaled”  denotes  timings  for  runs  in  which 
problem  size  is  increased  in  direct  proportion  to  number 
of  processors.  The  left  pane  of  Figure  4 shows  fixed 
system  size  speedup  for  CP2K  on  various  HPC  systems 
for  simulation  cells  of  different  shapes  and  sizes.  Each 
simulation  cell  is  composed  of  a series  of  unit  cells 
replicated  in  the  three  Cartesian  directions.  Each  unit  cell 
is  composed  of  8 nitrogen  atoms  in  the  cg-N  crystal 
structure.  The  shape  and  size  of  the  simulation  cell  is 
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denoted  by  h x k x 1,  where  h,  k and  I are  integers. 
Sapphire  (Cray  T3E)  produces  the  best  fixed-size  speedup 
with  a large  number  of  processes  for  near-cubic  systems; 
JVN  (Linux  Networx  Cluster)  exhibits  much  worse 
scaling  behavior.  However,  for  simulation  cells  that  are 
extremely  large  in  one  dimension  relative  to  the  other  two 
(such  as  those  used  in  our  shock  wave  simulations),  JVN 
demonstrates  greatly  improved  scaling.  On  JVN  we  see 
scaling  of  about  76%  of  ideal  when  doubling  the  number 
of  processors  for  our  simulations  up  to  256  processors, 
and  less  than  64  processes  hits  a memory  swapping 
bottleneck.  The  right  pane  of  Figure  4 shows  the  scaled 
speedup  for  DFT  calculations  using  CP2K  on  various 
HPC  platforms.  DFT  is  an  ~0(N3)  problem;  thus,  we 
expect  the  optimal  scaled  speedup  to  appear  ~0(N2).  As 
is  evident  from  the  right  pane  of  Figure  4,  we  see  Jaws 
(Dell  Woodcrest  Cluster)  outperforming  other  platforms. 
As  in  the  fixed  size  timings,  JVN  is  the  worst  performer. 
The  relative  performance  between  platforms  is  expected 
from  the  speeds  and  types  of  the  processors.  The 
relatively  poor  performance  of  JVN  in  both  sets  of  these 
timings  can  be  attributed  somewhat  to  the  age  of  the 
machine.  Consequently,  we  expect  the  performance  of 
the  Army  Research  Laboratory’s  (ARL’s)  new  MJM 
(Linux  Networx  Woodcrest  cluster)  to  exceed  that  of  Jaws 
and  the  other  systems. 

6.  Summary  and  Conclusions 

The  results  that  have  been  generated  show  no 
anomalous  behavior  and  we  project  no  unforeseen 
difficulties  in  achieving  the  scientific  goals  of  this  project. 
We  expect  to  conclude  the  cg-N  shock  calculations  by 
the  end  of  July  2007  and  are  currently  designing  the 
simulation  cells  that  will  be  used  in  shock  simulations  of 
the  conventional  carbon-hydrogen-nitrogen-oxygen 
energetic  material,  nitromethane.  At  this  time,  we  have 
simply  relaxed  a large  simulation  cell  of  nitromethane  to 
establish  the  template  that  will  be  used  to  generate  the 
filament  that  will  ultimately  be  subjected  to  shock.  This 
template  is  composed  of  3,584  atoms  arranged  in  the  low- 
temperature,  ambient  pressure  state.  NVT-MD 
simulations  of  this  system  have  used  about  50K.  hours  on 
JVN. 
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Figure  2.  Snapshots  of  cg-N  undergoing  a phase 
transition  at  t = 0 (left-most  frame),  t = 350  fs  (middle 
frame)  and  t = 700  fs  (final  frame)  at  T = 4,000K,  P = 
50GPa 
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Figure  3.  Snapshots  of  left  half  of  the  equilibrated 
crystal  at  the  beginning  of  the  shock  simulation 
(upper)  and  at  1.0  ps  after  shock  compression  (lower) 


Figure  4.  Scaling  and  Performance  of  CP2K.  The  left 
pane  shows  the  fixed  system  size  speedup  and  the 
right  pane  shows  the  scaled  speedup  for  various 
systems  and  sizes. 
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